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1. Introduction 

Precise calculations of the matrix elements of iso-scalar and iso-vector bilinear quark opera¬ 
tors within nucleon states are needed to probe many exciting areas of the Standard Model (SM) 
and its extensions. In Ref. [1], we showed that new scalar and tensor interactions at the TeV scale 
could give rise to corrections at the 10 3 level in precision measurements of the helicity flip parts 
of the decay distribution of (ultra)cold neutrons (UCN). This sensitivity is reachable in experi¬ 
ments currently under construction and being planned. Even if these experiments see a signal, to 
constrain the allowed parameter space of beyond the SM (BSM) models, however, requires that 
matrix elements of isovector scalar and tensor bilinear quark operators are known to 10-20% ac¬ 
curacy. Similarly, in Ref. [2], we showed that to probe novel sources of CP violation in neutron 
electric dipole moment (nEDM), a combination of matrix elements of the iso-scalar and iso-vector 
tensor operators are needed to estimate the contribution of the quark EDM to the nEDM. Lattice 
calculations are well poised to provide these estimates with the desired precision. 

In these proceedings, we present a detailed analysis of statistical and systematic errors in 
such calculations using 9 ensembles of 2+1+1 flavor HISQ lattices generated by the MILC col¬ 
laboration [3]. The matrix elements are calculated using clover valence quarks. We examine the 
following sources of systematic errors - contribution of excited states, estimates of renormalization 
constants, finite volume and lattice discretization effects and dependence on quark mass. Three of 
these sources, statistics, contribution of excited states, and renormalization constants, affect the 
precision with which estimates from an individual ensemble are extracted. The other three, finite 
volume and lattice discretization effects and dependence on quark mass, require fits and extrapola¬ 
tions based on all the points. We examine the two classes of uncertainties separately. 

2. Statistics 

The MILC Collaboration [3] has generated ensembles of roughly 5500 trajectories of 2+1+1- 
flavor HISQ lattices at three values of light quark masses corresponding to M n ~ 310, 220, 130 
MeV at a = 0.12, 0.09 and 0.06 fm. We analyze configurations separated by 4-6 trajectories of the 
hybrid Monte Carlo evolution and discard the initial 300-500 trajectories for thermalization. The 
status of our analyses using these ensembles are summarized in Table 1. To increase the statistics, 
each configuration is analyzed using gaussian smeared sources in multiple locations displaced both 
in time and space directions to reduce correlations. 

We performed the following statistical tests. The data for a given ensemble are divided into 
bins (by source point and configurations) and the Kolmogorov-Smirnov test is performed on quan¬ 
tities that have reasonable estimates configuration by configurations (iso-vector vector charge, 
value of 2-point function at a given time separation). While, this pairwise test showed that the 
sub-samples are consistent with being drawn from the same distribution, the mean values of ob¬ 
servables fluctuated by up to 3a. This variation is much larger than expected based on our bin size 
of over 1000 measurements. We do not find long tails in the distributions for any of the samples 
that could explain the fluctuation, but do see a variation in the sample distribution. One possible 
explanation is that the ensembles have not covered enough phase space and errors are consequently 
underestimated. Our overall conclusion is that a few thousand, and possibly 0(10,000) for the 
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Label 

L 3 xT 

M„ MeV 

(M n L) 

Acfgs 

N Measurements 

^sep 

al2m310 

24 3 x 64 

305.3(4) 

4.54 

1013 

8104 

8, 9, 10, 11, 12 

al2m220S 

24 3 x 64 

218.1(4) 

3.22 

1000 

24K (12K) 

10(8, 12) 

al2m220 

32 3 x 64 

216.9(2) 

4.3 

958 

7664 

8, 10, 12 

al2m220L 

40 3 x 64 

217.0(2) 

5.36 

1010 

8080 

10 

a09m310 

32 3 x 96 

312.7(6) 

4.5 

881 

7058 

10, 12, 14 

a09m220 

48 3 x 96 

220.3(2) 

4.71 

890 

7120 

10,12,14 

a09ml30 

64 3 x 96 

128.2(1) 

3.66 

883 

4824 

10, 12, 14 

a06m310 

48 3 x 144 

319.3(5) 

4.51 

865 

3460 

16,20 

a06m220 

64 3 x 144 

229.2(4) 

4.25 

430 

1320 

16, 20, 22, 24 


Table 1: Description of the nine ensembles at a = 0.12, 0.09, 0.06 fm used in this study. 


scalar charge, independent configurations with 0(32) measurements on each are needed to obtain 
estimates with < 2% errors. 

3. Excited-State Contamination 

Our current data show significant excited state contamination in both 2-point and 3-point func¬ 
tions. The goal is to extract all observables (charges, charge radii, form factors) by calculating 
matrix elements between ground-state nucleons. We address this by using smeared operators tuned 
to increase coupling to the ground state and suppress radially excited states. Second, as discussed 
in [4], we partially remove the remaining contamination by including one excited state in the anal¬ 
ysis. Higher states are not included because with current statistics we are not able to resolve them, 
especially in the 3-point functions. 

Denoting the first excited state mass by M\ and coupling to our operator by sf\ , the three-point 
function with source at f,- = 0, operator insertion at t = t and sink at tf = t sep can be written as 

\^\ 2 {0\O r \0)e- M ^f-^ + |M | 2 (11O r 11) e~ Ml 
+ (01 O r 11) e~ M ° g _Ml ( f / _f) + 

+ ,<M(l|Or|0)e" Ml(( " f 'V Mo ^~ f) . (3.1) 

The masses and amplitudes Mo, M\, My, and si\ are obtained from fits to the two-point functions. 
The desired matrix element (0|Or|0) is then obtained by isolating (0|Or|l) and (I |0i j l). This 
requires doing calculations with multiple values of t and t sep . Using the sequential source method, 
we carry out operator insertion at all values of t between the source and sink timeslices. The values 
of t sep investigated are listed in Table 1. A nonlinear least-square fitter is then used to extract 
(0|Or|0) by fitting the data for all f sep simultaneously using Eq. (3.1). 

The a = 0.12 data for all three charges do not exhibit significant trends with respect to t sep [4], 
Simultaneous fits to f sep = 8, 10 and 12 data are consistent with a fit to just the t sep = 10 data. 
Consequently, in [4] we had concluded that f sep > 1.2 fm is needed to control excited state con¬ 
tamination. The a = 0.09 and 0.06 data are much cleaner and show an increase in with ? sep as 
illustrated in Fig. 1. On the other hand gj continues to show very little senstivity to f sep . The errors 
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Figure 1: Fit using Eq.(3.1) to extract unrenormalized gA.s.r from the a09m220 ensemble data. The black 
line and the grey error band are the result of the simultaneous fit to the f sep = 10, 12 and 14 data. The 3 colored 
lines are obtained from the simultaneous fit evaluated for each f sep . To reduce excited-state contamination, 
the grey points on either side that are close to the source/sink timeslice, are not included in the fit. 



Figure 2: Study of finite volume effects in unrenormalized gA,s,r using fits to Eq.(3.1) on the f sep = 10 data 
from the three al2m220 ensembles. Rest is same as in Fig. 1. 


in gs are too large to draw conclusions. Overall, trends in gA data at weaker couplings imply that 
4ep > 1.5fm is needed to establish control over excited state contamination. 

Our fits are biased by the smallest f sep data because the statistics are the same for all f sep , while 
errors increase with f sep . For example, on the a = 0.12fm ensembles, the statistical errors increase 
by about 40% with each unit increase in f sep . This estimate scales with a, i.e., on the a = 0.06fm 
ensembles, the same fractional increase takes place every 2 units. Similarly, the errors increase 
by about 20% on lowering the light (u and d) quark masses by a factor of two, i.e., going from 
M n = 310 to 220 MeV ensembles. As an illustration, consider fits to a = 0.09fm ensembles with 
fsep = 10, 12, 14 shown in Fig. 1. The f sep = 10 data make the largest contribution to the extraction 
of (0|Or|0) and (0|(9|j I). The change between ? sep = 10 and 12 contributes to fixing (1 |Op| 1). For 
f se p = 14 data to contribute at the same level, its statistics should be 3-4 times that of f sep = 10 data. 

4. Finite Volume Effects 

The results of our finite volume study using the al2m220 ensembles with volumes 24 3 , 32 3 and 
40 3 (corresponding M K L = 3.22, 4.3 and 5.36) are shown in Fig. 2. The gA data show significant 
increase with volume, while the gj data show saturation between the M n L = 4.3 and 5.36 data. 
The gs data are again too noisy. Our conclusion is that lattices with M n L > 5 are needed to control 
finite volume effects unless we can reliably model extrapolation in lattice volume. 
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5. Non-perturbative Renormalization 

We use the RI-sMOM scheme to calculate the renormalization constants of the bilinear quark 
operators [5], Details of the method and our implementation were presented in [4]. The most im¬ 
portant issue, especially when using smeared lattices, was demonstrating the presence of a window 
in momentum q, A qcd <C g <Cc /a with e an a priori unknown number of 0 ( 1 ), where lattice arti¬ 
facts are expected to be small. Sufficiently close to the continuum limit where perturbation theory 
works, a test of whether such a window exists is that Zg and Zy in the RI-sMOM (or any lattice) 
scheme should show a q 2 dependence given by the anomalous dimensions of these operators along 
with a weaker dependence from the running of cc s , while Z A should only show the latter. Con¬ 
verting estimates obtained from within such a window in q 2 in the RI-sMOM scheme to the MS 
scheme used in pehnomenology and run to some fixed scale, say p = 2 GeV, should give estimates 
independent of q 1 . Based on our analyses at the three lattice spacings (see [4] for details) our con¬ 
clusions are: there is evidence of such behavior, i.e., a window, in Z,\ and Zy, but not in Z$ even 
on a = 0.06 fm ensembles. Lacking a convincing demonstration of a window, we have defined 
a procedure that will extrapolate to the right continuum limit [4] and have been conservative in 
estimating errors, however, we recommend a further study at various a, in particular for Z 5 . 

6. Combined fits in lattice volume, spacing and quark mass 

Having discussed the first class of uncertainties that affect individual data points, we now 
discuss extrapolations in lattice spacing and volume, and the quark mass. It is very hard to generate 
dynamical lattices with fixed quark masses (fixed Mf) and lattice volumes (fixed M n L ) at multiple a 
in order to take the continuum limit along a line of constant physics. Similarly, it is not easy to hold 
the lattice volume and a constant and vary the quark mass to study the chiral behavior. Our best 
option is to do a combined fit in a , M n and M n L to obtain physical estimates. The second challenge 
is choosing the extrapolation ansatz in each of these three variables — we have to compromise 
between the number of free parameters included and the number and quality of data points. In 
Fig. 3 we show such a fit keeping only the leading order terms in each of the three variables, 

g{a,M Kl M n L) = g physical + aa + /3 m\ + ye ~ M * L . ( 6 .1) 

In Fig. 3, note that the errors in individual points vary significantly. As a result, with 9 data points, 
this ansatz with 4 free parameters is the most extensive we can explore. 

Fig. 3 summarizes the trends mentioned before. Removing excited state contributions and 
doing finite volume and chiral extrapolations all increase gA towards the experimental value. Data 
for gr show almost no dependence on cl M k or M n L and give gj = 1.06(0.06). We consider this 
estimate reliable. Statistical errors are too large to draw conclusions about g$. 

Prognosis for the future'. Based on current analyses, we conclude that with current ensembles 
and (9(10,000) measurements, we can extract gj with about 2% errors on each point, and ~ 5% 
uncertainty in the extrapolated value. To extract g A with similar precision will require (9(2000) 
configurations at three values of a < O.lfm, M n L > 5 and (9(24) measurements on each configura¬ 
tion to get statistically significant data for 1.2 < t sep < 1, 6 fm. Estimates of gs with similar precision 
will require an order of magnitude more measurements. 
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Figure 3: Results of combined fits in a, M n and M n L to obtain the physical results for the renormalized 
charges using the ansatz given in Eq. (6.1). Iso-vector charges are labeled as g" d . 


7. BSM contributions to Neutron Electric Dipole Moment 


Lattice calculations of matrix elements of effective quark EDM and chromo EDM operators 
within a neutron state to proble BSM theories were initiated in [2]. The simpler is the quark EDM 
which is an extension of gT but matrix elements of both isovector and isoscalar tensor operators 
are need. One, therefore, has to evaluate and control the signal in the disconnected diagrams [6]. 
We also analyze operator mixing and renormalization in 1-loop perturbation theory. For brevity, 
operators that vanish by the equations of motion are included by introducing the field combinations: 

Ve = 0 iD M Yg ~ m)w , = dn ~ igA^T* - ie w A^ (7.1) 

y/E = — t/f +m) , ^ = d ^ + igA‘^ T a + ieyA^P . (7.2) 

In terms of these fields, the operators we study are given in Table. 2. The pattern of mixing of the 
dimension 5 operators under renormalization that needs to be calculated is given in Table 3. Papers 
containing 1-loop results for the mixing and a first estimate of the quark EDM are being prepared. 
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(?( 3 ) = iP = yiy, 

of = C='{\j/(a^ v Y 5 + Y 5 a^ v )G^ v \i/ 

Of = id 2 P 


of = E = f ^ r (a^ v 75 + 7 5 a /JV )i^ v v/ 

of=GG=^ va PG^ ap 

of = m FF 

of = d-A = dniyrfTfcY) 

Of = m GG 

of = imP = miff iff, 

of = md A 

of =FF = \e^F^F ap 

Of = m 2 iP 

of = IFee. = i Ye 75 Ye 


Of = d-A E = d,j [\j/ E Y f "Y5V+ VY^YsVe] 

Of =A d = y/Ys^VE ~ Ye $ YsY 


0 ( n =A Mr) = ie (y4. (y) Y5VE - Ve/^Ysv) 


Table 2: Operators of dimension 3, 4 and 5 needed in the nEDM calculation. 



C 

d 2 P 

E 

mFF 

mGG 

md ■A 

m 2 P 

Pee 

d - Ae 

A d 

Aam 

C 

Z c 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

d 2 P 

0 

Zp 

0 

0 

0 

0 

0 

0 

0 

0 

0 

E 

0 

0 

Zt 

0 

0 

0 

0 

0 

0 

0 

0 

mFF 

0 

0 

0 

Z m l Z F F 

0 

0 

0 

0 

0 

0 

0 

mGG 

0 

0 

0 

0 

Z m 1 Z G Q 

X 

0 

0 

0 

0 

0 

md A 

0 

0 

0 

0 

0 

Z m 1 Z-)A 

0 

0 

0 

0 

0 

m 2 P 

0 

0 

0 

0 

0 

0 

7 1 

0 

0 

0 

0 

Pee 

0 

0 

0 

0 

0 

0 

0 

X 

X 

X 

0 

d -A e 

0 

0 

0 

0 

0 

0 

0 

0 

X 

0 

0 

A d 

0 

0 

0 

0 

0 

0 

0 

X 

X 

X 

0 

A Air) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

X 


Table 3: Mixing due to QCD of the dimension-5 operators. Non-zero entries need to be determined. 
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